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We devise a non-perturbative method, called Parametric Perturbation Theory (PPT), which is 
alternative to the ordinary perturbation theory. The method relies on a principle of simplicity for 
the observable solutions, which are constrained to be linear in a certain (unphysical) parameter. 
The perturbative expansion is carried out in this parameter and not in the physical coupling (as in 
ordinary perturbation theory). We provide a number of nontrivial examples, where our method is 
capable to resum the divergent perturbative series, extract the leading asymptotic (strong coupling) 
behavior and predict with high accuracy the coefficients of the perturbative series. In the case of 
a zero dimensional field theory we prove that PPT can be used to provide the imaginary part of 
f^*) . the solution, when the problem is analytically continued to negative couplings. In the case of a 

f*"^ ' 4 lattice model 1 + 1 and of elastic theory we have shown that the observables resummed with 

£N) | PPT display a branch point at a finite value of the coupling, signaling the transition from a stable 

to a metastable state. We have also applied the method to the prediction of the virial coefficients 
for a hard sphere gas in two and three dimensions; in this example we have also found that the 
solution resummed with PPT has a singularity at finite density. Predictions for the unknown virial 
coefficients are made. 
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For the majority of the problems in Physics no exact analytical solution is known and it is a common procedure 
to resort to Perturbation Theory (PT) [J H, H, 3 ■ The fundamental idea behind perturbation theory is that, if a 
given problem is solvable for a particular value of a parameter, then one can obtain analytical approximations to the 
solution in a close neighborhood of this value, by Taylor expanding in that parameter. Calling g this parameter and 
^3 ' 9o the value for which an exact analytical solution is known, then the results obtained with perturbation theory will 
be, to a given order, polynomials in (g — go). The perturbative series obtained by considering all the terms of the 
expansion will in general have a finite radius of convergence, r, which in some cases could even be zero and therefore 
the series would be divergent for all g ^ 0. Actually divergent series are usually expected from the application of 
■ PT to quantum field theory, as first observed by Dyson in the case of Quantum Electrodynamics (QED)j5|]. Another 
\ well-known example of divergent series is given by the quantum anharmonic oscillator, whose perturbative coefficients 
for the energy of the ground state have been calculated by Bender and Wu in Q and proved to have a factorial growth. 

Although the pertubative series provide in many cases the only systematic approach to the solution of a problem, 
they are not always useful, since they are confined to a restricted region for the physical parameters, \g — go\ < r. 
Outside this region, the physics becomes nonperturbative and cannot be described directly in terms of the original 
series. For quite a long time physicists have been interested into finding a bridge from the perturbative to the non- 
perturbative region. Several methods have been developed which allow to extract the non-perturbative behavior from 
the perturbative series: among such methods we would like to mention the Borel and Pade approximantsjl], 0] and 
nonlinear transformations 0]. 

Methods which are alternative to PT should retain on one hand the ability to provide a sistematic analytical 
approximation to a given problem, and on the other hand they should remain valid even in the nonperturbative 
region, never leading to divergent series. Over the years new ideas have allowed to devise methods which comply with 
these requirements. The Linear Delta Expansion (LDE) Q and the Variational Perturbation Theory (VPT) [l(| are 
two examples of non-perturbative methods. Roughly speaking these methods work by introducing in a problem an 
artificial parameter and turn the original problem into a new one with a modified perturbation. The optimization 
of the "perturbative" results to a given order with respect to the artificial parameter is usually obtained through 
the Principle of Minimal Sensitivity (PMS) [ll[ and leads to expressions which are non-polynomials in the physical 
parameters and therefore non-perturbative. 
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In this paper we will explore a new a path, which is also described in shorter letter: the method that we have 
devised, which we have called Parametric Perturbation Theory (PPT) method, is based on few simple ideas. The 
first one, which we will refer to as Principle of Absolute Simplicity (PAS), is that we do not want to calculate the 
observable (energy, frequency, etc.) directly as a polynomial in the physical coupling g, as done in PT, but that 
this observable should have the simplest possible form (linear) in a given unphysical parameter g; the second idea is 
that the perturbation theory must be carried out in g and that the functional relation g = g(g) must comply with 
the Principle of Absolute Simplicity to the order to which the calculation is done. This will allow to determine the 
relation between g and p and in turn to obtain the observable as a parametric function of g. 

The paper is organized as follows: in Section [II] we develop the Parametric Perturbation Theory for a problem of 
nonlinear oscillations in classical mechanics and show that at finite order it provides extremely accurate approximations 
for the frequencies and for the solutions of the problem; in Section HlTl we show that the PPT approach can be applied 
directly on the perturbative series and discuss the performance of our method in the case of non trivial examples, 
with divergent perturbative series. Finally, in Section HVl we draw our conclusions. 



II. THE METHOD 



Consider a model which depends on a parameter g, and which is solvable when g = 0. The application of PT to this 
problem to a finite order yields a polynomial in g. Calling r the radius of convergence of the perturbative series, the 
direct use of PT must be restricted to \g\ < r, as previously discussed. However, the misbehavior of the perturbative 
series for a physical observable O is the result of having expanded in a parameter, <?, which is not optimal. If one 
knew the exact solution to the problem, i.e. O = /(<?), then this solution could be considered as a polynomial of order 
one in the variable g = f(g). Although this observation by itself cannot be used as a constructive principle, we may 
adopt the philosophy that the perturbative series for the observable can be simpler and convergent in all the domain, 
if it is cast in terms of a suitable parameter g. 

Only if such parameter, by luck or ability, turns out to be the g = f(g) discussed above, the exact solution is 
obtained. The goal, therefore, is to progressively build this parameter g to yield an expression for O as simple as 
possible. In this framework the perturbative expansion is carried out in g and all the physical quantities in the 
problem are expressed as functions of g. In particular we have now that g = g(g). While the ordinary perturbation 
theory works by calculating the contributions to higher orders in g 7 each term of higher order refining the result to 
lower order, the approach approach is the opposite: we carry out a perturbative calculation in g, and then determine 
order by order the form of g = g(g) so that the observable O(g) can be a order one polynomial in g. This is in essence 
the Principle of Absolute Simplicity. 

Having given the general ideas of the method we proceed to examine its implementation in a concrete problem. We 
consider the classical nonlinear oscillations of a point mass described by the equation (Duffing equation) 

^-+ x (t) = -gx 3 (t) . (1) 

The Lindstedt-Poincare method can be used to obtain a perturbative expansion of the squared frequency f2 2 of the 
oscillations in powers of g [H, 0, Q- The method works by defining an absolute time scale, independent of g and by 
then fixing the coefficients of the expansion of fl 2 so that the secular terms in the expansion are eliminated at each 
order. Working through order (gA 2 ) 5 one finds 

o2 1 3gA 2 3g 2 A A 9g 3 A 6 1779 ff 4 A 8 ^ r . /l2 . 51 

n 2 Kii + — - v— - \-o \{qA 2 y\ , (2) 

4 128 512 131072 W 

A being the amplitude of oscillations 1 . 

We now proceed to implement our method. The first step is to define a functional relation between g and the 
perturbative parameter of the expansion, g. For example, we choose 

9(0) = e 1+ gr c " g " ' ^ 



1 The radius of convergence of the perturbative series in this case is g = 1/A 2 . 
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where the coefficients c n and d n are unknown constants to be later determined. The parameter N is related to the 
order to which the calculation is performed, since the number of coefficients c, and d^ needs to match the number of 
conditions available at a given order. 

The reader may wonder the reason of the particular choice made in eq. Q : the functional relation between g and 
g can be more general than eq. ([3]), although it is not completely arbitrary since it must reproduce all the terms in 
the perturbative expansion when the relation between g and g is inverted. The choice made here takes into account 
this fact and also the leading asymptotic behavior of the frequency lim^oo tt 2 oc g? 

We now follow the standard procedure of the LP method and introduce an absolute time r = fl(g) t. Applying the 
PAS we choose this relation to be 

n 2 (p) = a 1 +a 2 g , (4) 
where ai,2 are coefficients to be determined. We also expand the solution as 

oo / \ n 

y(T) = Yv n (T) . (5) 



n=0 



<J 



and the solution 



The reader familiar with the LP method will recognize the profoundly different character of the present approach: 
in the LP method the observable, i.e. f2 2 , is expressed as a series in powers of g, whose coefficients are determined by 
the condition that secular terms are eliminated at each order; here we impose that Q 2 has the simplest possible form 
when expressed in terms of g, and we let g{g) to contain arbitrary powers of g. In detail, we transform the original 
equation into the new equation 

n 2 (g)^+y(T) = -g(g)y 3 (T) . (6) 

For sake of simplicity we will limit ourselves to work through order g 3 and solve the differential equations resulting 
at each order in g. The elimination of the secular term to order one yields ol\ = 3A 2 /4, as in standard LP method 
(cto = 1). The solutions to order and 1 are Ho(t) = A cost and y\{r) = 32 ( c f „rf 1 ) ( — cost + cos3r). 

The elimination of the secular term to second order provides the condition 

(7) 

y 2 (r) = ^23 A - cos t - (2AA - 5^ cos 3r . (8) 

Finally, one can determine c\ to cancel the secular term to third order, c\ = || A 2 , and correspondingly d\ = j^-A 2 . 
The solution to third order thus reads 

?/3 (t) = 5 A cos r — AA cos 3r — 2 A cos 5r + A cos 7r . (9) 

To order g 3 we therefore find 

_ 32 + 23gA 2 

9 W 9 32 + 22gA 2 ' (10) 
which can be inverted and used to express the frequency directly in terms of g: 

Q 1 1 33 3 , 

n 2 « 1 + - Q A 2 = - + -gA 2 + -^(121.9^ + 384)^+ 256 . (11) 

Notice that this expression is non-perturbative in gA 2 and that it provides a maximum error E = 
lim fl ^2_ >00 (l — Q 2 /n 2 xact ) fa 5.27 x 10~ 4 . This error is much smaller that then the one obtained to third or- 
der in [l2| using the Linear Delta Expansion, i.e. &lde = 699 " 4 96gy^'+i28 +128 ( ul wmcn ca se one has Slzjb = 
Hindoo (1 - n 2 LDE /n 2 exact ) « 1.36 x 10- 2 ). 



2 We are assuming that the denonimator of eq. does not have zeroes for g > and that therefore g = 00 is reached for q = 00. 



4 




1U()() 



FIG. 1: (color online) Left panel: Error over the square frequency, E = 1 — Qapprox/^lxactj calculated to order 3,5 and 7 (solid, 
dashed and dot-dashed curves); Right panel: Error over the first two Fourier coefficients calculated to order 3 using PPT and 
the LDE approach of [f3} |. 



At the same time PPT provides highly accurate estimates for the Fourier coefficients of the solution, c n . In the 
right panel of fig[T]we have plotted the error defined as H n = \ c c ^PV rox j c exact _ i| a s function of g and compared our 
results with the already excellent results of EH . 

At this point the reader should recognize that, if one is interested only in frequency and not in the solution, it is 
possible to determine the coefficients c„ and d n to a given order with a minimal effort directly from the coefficients of 
the pcrturbative series. The procedure consists of first substituting g = g(g) inside the perturbative series, and then 
expanding around g = 0: the unknown coefficients c„ and d n are then used to suppress the nonlinear behavior in g 
inside f2 2 . Following this simple procedure we have produced the result to order 5 and 7 in the left panel of FigQ] 



III. RESUMMATION OF PERTURBATIVE SERIES 



In many cases the perturbative results for a given problem are known only to a finite order since the calculation 
of each higher order involves increasing technical difficulties. For example, the calculation of observables in Quantum 
Field Theory at a given order in PT requires to take into account a number of Feynman diagrams which rapidly grows 
with the order of the perturbation, the calculation of the higher order (multiloop) diagrams being more and more 
challenging. 

For this reason it is desirable to have a procedure which, using only a finite number of perturbative coefficients, 
may extract the essential physical behavior of the solution and possibly predict a number of unknown perturbative 
coefficients. We will show now that this result can be efficiently achieved using our method. 



A. Anharmonic oscillator 



Consider the quantum anharmonic oscillator 

£ = y + y+^ 4 - (12) 

The series obtained with perturbation theory for this problem is divergent and its coefficients behave as 

&„« (-l) n+i y67^r(n + 1/2)3" (13) 

for n — > oo as shown in Q. These coefficients can be also obtained exactly, using the recursion relations given by 
Bender and Wu in J6[. The resummation of this perturbative series has been considered by several authors, using 
different techniques Q, i, EI EE EE EE El EE HE lU, H3, H| . 

We will follow the philosophy of our method and define a functional relation between the physical coupling g and 
the unphysical parameter g. In principle this relation can be expressed by mean of an arbitrary function, but it can 
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TABLE I: Comparison of the exact perturbative coefficients for the anharmonic oscillator with the approximate coefficients 
obtained with PPT using N = 3. 
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also use information coming from the strong coupling regime, where the energy goes like Ea oc g 1 / 3 as g — > oo. For 
example we can choose 



gio) = q 



(14) 



which has the correct asymptotic behavior for g ~ > oo, provided that the denonimator does not vanish for g > 0. The 
unknown coefficients c n and d n in this expression will be determined so that the ground state energy is linear in the 
unphysical parameter, as required by the PAS, i.e. 

E =bo + b ig . (15) 

Choosing A = 3, corresponding to use only the first 6 perturbative coefficients, we can fully determine the coefficients 
c„ and d n : 

3111725471 292194444505 1136953355311 



109345364 1749525824 J 6998103296 

730092771 _ 215945995035 

1 _ 27336341 ' 2 ~ 1749525824 ' 
Working to this order it is possible to extract the leading coefficient of the strong coupling series 

b + b ig ( 215945995035 \ 2/3 . . 

a = lim . , = 3 —————— — — w 0.624458 , (16) 

e ^oo g (g) V 2273906710622 / 

which should be compared with the fairly precise results of [J H3, H3, a ° = 0.66798625915577710827096. In the 
opposite limit, g — > 0, the energy calculated with the PPT can be cast in terms of g after inverting 



En 



3.9 21g 2 
T 8~ 
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0\g 



hi 



(17) 



which is correct up to order g e . In Table U we compare the exact perturbative coefficients going from the order g 7 
with the approximate ones predicted by the PPT using N = 3. The last column displays the error 



to order g 
E„ = 100 x 



/vexact 
On I V n 



1 



The reader could argue that the quality of the results that we have obtained is due to having taken into account 
the exact asymptotic behavior of Eo for g — > oo. We will now show that our method allows one to optimize the 
asymptotic behavior of the approximate solution, even in the case where the exact behavior is unknown. We consider 
the approximations corresponding to 



9{q) = Q 



l + En=lCng" 



1 2 



(18) 
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TABLE II: Comparison of the exact perturbative coefficients for the anharmonic oscillator with the approximate coefficients 
obtained with PPT using N = 3. 
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fixing N u + iVrf = 5 as in the previous case. We can compare the first coefficient predicted by our method, 67, using 
the different sets. These coefficients are displayed in the first row of Table [TT1 the second row displays the error (in %) 
with respect to the exact coefficient, shown in Table HI The set [3,2] previously considered provides the lowest error 
and therefore selects the correct asymptotic behavior. 

This example shows that, even in the unfortunate case where the asymptotic behavior of the energy is unknown, 
it is possible to extract some information on the strong coupling regime directly from the perturbative series, using a 
limited number of perturbative coefficients. 

In Fig[2] we have plotted the exact energy (numerical) as a function of g (squares) and we have compared it with 
the results of PPT applied to 3 different orders, all reproducing the exact asymptotic behavior. Our results approach 
the numerical result as N is increased. 



B. A PT symmetric hamiltonian 



The complex hamiltonian 



H 



igx 



(19) 



has been the first example of a PT symmetric hamiltonian which has a completely real spectrum to be discovered. 
Bender and Dunne have studied in 25[ the large order perturbative expansion of the ground state energy of this 
hamiltonian finding that the coefficients of this series grow as 



b n « (-1) 



n+l ou 

(2tt) 3 / 2 



r(n+l/2) [l-0(l/7i)] 



(20) 



Table I of |25[ contains the first 20 coefficients of the perturbative series. In a recent paper Bender and Weniger [261 ] 
have provided numerical evidence that the perturbative series for this PT symmetric hamiltonian is Sticltjcs, using 
the first 193 nonzero coefficients. 

We will here use this model to obtain a further test of PPT. We assume the functional relation 



;(«) 



1 + En=l C nQ 



(21) 



where the difference N — M constrains the asymptotic behavior, which is not known exactly in this case. Notice 
that the square root in the definition of g is a consequence of the fact that the perturbative series contains only even 
powers of g [25|. In FigUwe have compared the exact numerical results of the last column of Table III of [25| with 
the calculation obtained with PPT using three different sets of (N,M), which correspond to the same number of 
conditions 3 . Our results suggest that the asymptotic behavior of the energy is approximately E$ cx ^fg for j->oo. 

C. Zero dimensional 4> 4 theory 



Integrals of the form 



r+00 

E{g) = / e- x2 - gx dx 
Jo 



(22) 



Of course, when numerical results are not available, one can resort to the same approach followed for the anharmonic oscillator, selecting 
the optimal asymptotic behavior among those available. 
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FIG. 2: (color online) Ground state energy of the anharmonic oscillator as a function of g. The squares are numerical results, 
the curves correspond to the results obtained with PPT to different orders. 




FIG. 3: (color online) Ground state energy of the anharmonic oscillator as a function of g. The squares are numerical results, 
the curves correspond to the results obtained with PPT to different orders. 



can be used as a model of a 4 in zero dimensions [13, [H, H^. As for the case of higher dimensional field theories, the 
perturbative scries for this model is divergent. The authors of [28l. |29| have proved that the Linear Delta Expansion 
(LDE) is able to deal with this problem and that it provides results which rapidly converge to the exact value. 
The integral in eq. (|22p admits an exact analytical solution which is given by 



E{g) 



#1/4 



1 



(23) 



where Kiu(g) is the Bessel function of order 1/4. Notice that for negative values of g this expression acquires an 
imaginary part, signaling that the system becomes metastable. 

We will now analyze this problem with the help of PPT. We choose the functional relation 



g(e) 



i + E 



A' 



-I 5 



C n Q 



which allows one to obtain the correct asymptotic behavior, E(g) oc g~ 1//4 



(24) 



as g 



CO. 
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As usual the application of the PPT method requires that the coefficients c n and d n corresponding to a given N 
be determined by imposing the PAS, i.e. by asking that the observable E{g) be linear in the parameter g. Using 
three different sets, corresponding to JV = 1,2 and 3 we have observed that our results converge quickly to the exact 
analytical result (see the left panel in Figf?]). 

We will however move further and concentrate over the best set, N = 3. In this case the relation between g and g 
is given by 

p(2924.98£ 3 + 881.78 £ > 2 + 58.83 (? +l) 5 , N 

gig) ps r . (25) 

(-1737.20p 4 + 2243.93g 3 + 832.17(? 2 + 57.96g + l) 5 

Since go = —0.0593 is a zero of the denominator, lim e _> 0[l g(g) = oo: this result signals the presence of a branch point 
in the proximity of go (see the right panel of Fig|4]). We now consider the region g < 0, where the analytic continuation 
of the solution acquires an imaginary part. Using eq. (|25j) we find the numerical solutions of the equation g(g) = g, 
with g < 0. For example, corresponding to g = — 1 we find two pairs of complex conjugated roots accompanied by a 
single real root: 

qi = 0.20784008231963882+ 0.5489736369789899 i (26) 
02 = 0.20784008231963882 - 0.5489736369789899 i (27) 
g 3 = -0.24838105054544726 (28) 
g 4 = -0.2232185118896695 + 0.006357863317100693z (29) 
g 5 = -0.2232185118896695- 0.006357863317100693i . (30) 

It is important at this point to notice that obtaining a complex value for g has an immediate effect on the observable 
E(g), which acquires an imaginary part, ImE(g) = b\ Img. We can verify if one of these solutions corresponds to 
the exact solution of (|22"|) for g = —1: 



ImE(g) = -0.3767931291206198 (31) 

which should be compared to the imaginary parts calculated with the PPT using the numerical roots gi, i = 1, . . . , 5: 

ImE(g) PPT = -0.36488641384088155 (32) 

ImE(g) PPT = 0.36488641384088155 (33) 

ImE{g) PPT = (34) 

ImE(g) PPT = -0.004225882244972266 (35) 

ImE{g) PPT = 0.004225882244972266. (36) 

These results suggest that the first root corresponds to the analytic continuation of the solution for g > to 
negative values. On the other hand the real part of ReE(g) PPT = 0.7480818175977717 has the opposite sign of 
ReE(g) = —0.7603309714715291: this happens because our function is continous and therefore it is not possible to 
reproduce a discontinuity at g = 0. 

To test the conclusions that we have just reached we can plot the real and imaginary parts of E(g) as obtained 
from (|22p and compare them with the results provided by the PPT. In FigJ5]we show the results obtained with this 
comparison. We conclude that both the imaginary and real parts of E{g) (apart for a sign) are reproduced with good 
quality; clearly the exponential behavior of the exact solution for g — > 0" cannot be reproduced in this approach. 

Let us now briefly explore a different issue. If we push forward the analogy with quantum field theory, the application 
of the PPT to this problem can have a simple interpretation in terms of Fcynman diagrams: because applying the 
PAS we are expanding not in the coupling, but in a parameter g, there is a infinite number of vertices appearing at 
tree level, whose couplings contain the unknown constants c n and d n . The spirit of the PAS is then to perform a 
perturbative (in g) calculation in which all diagrams corresponding to orders higher than g cancel out by fixing the 
unknown constants and therefore yielding the observable in terms of just two diagrams, with zero and one vertex 
respectively. This argument is sketched in Fig|6l 

D. QED effective action 

We consider the QED effective action in the presence of a constant background magnetic field: 



e 2 B 2 f°° ds f , 1 



.s 



9 





FIG. 4: (color online) Left panel: Comparison between the exact integral for the zero dimensional (j) 4 theory and three different 
approximations obtained using the PPT; Right panel: comparison between the set [3, 4] and the exact integral. The approximate 
solution has a branch point close to g — 0. 
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FIG. 5: (color online) Real and imaginary parts of E(g) for the zero dimensional (f> 4 theory obtained using the PPT with N = 3. 
The results are compared with the exact expression of eq, (|22|l . 



FIG. 6: Feynman diagrams for the propagator in a <f) A theory. The diagrams with dashed line need to cancel fixing the unknown 
coupling contained in the bold vertex. The coulings coming from the expansion of g to order g 2 and higher are represented 
with a bold circle. 
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FIG. 7: (Color online) Comparison between the exact numerical result for E{g) = S(g)/(2 e ^ ) and the approximation obtained 
with the set [4, 4]. 



where B is the magnetic field strength, e the electron charge and m e the electron mass. Following [30] ] we introduce 
the effective coupling g = e 2 B 2 jrr? e and obtain the divergent perturbative series 



P 2 R 2 00 
S = -2 6 —g Y^bn9 n 

n=0 



with 



(-1) 



ra+1 



4" IS 



2n+4\ 



(2n + 4)(2n + 3)(2n + 2) 



(38) 



(39) 



&2n+4 being Bernoulli numbers. 

In Fig[7]we show the comparison between the exact numerical result for E(g) = S(g)/(2 e ) and the approximation 
obtained with the set [4,4]: although we have used only ten perturbative coefficients, the resummation provides a 
quite precise approximation over a large range of the coupling in an analytical form. 



E. One plaquette integral 



In [31j the weak coupling expansion for a one-plaquette SU(2) lattice gauge theory was discussed. The partition 
function in this case is given by 



Z{fl) = - / Vl - U 2 e -/3(1~«) du 
K J-l 

and can be calculated exactly in terms of the modified Bcsscl function I\ : 

Z(f3) = 2e"^ 



(40) 



(41) 



This expression admits a convergent series expansion around /3 = (strong coupling expansion), but provides a 
divergent series when expanded in the opposite regime, f3 — > oo, (weak coupling expansion) [3 lj: 



(42) 



The terms of this series grow like l\/2 l and the sign does not oscillate. 
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FIG. 8: (color online) P versus j3 for SU(2) on one plaquette. The solid line is the exact result; the dashed line corresponds to 
the Weak Coupling Expansion; the dotted line corresponds to the Strong Coupling Expansion. 



We will now apply our method to this model, considering both regimes and using as usual the functional relation 



g(o) = q 

to be determined independently in the two regimes. 



1 + TZ=l C nQ' 



(43) 



Weak coupling expansion 



In this case we identify j3 = 1/g and use a set corresponding to N u = 2 and Nd = 3, obtaining the solution 

1 _ g(l.816g 2 - 3.428g+l) 

~ 9 ^ Q ' ~ 0.154g 3 + 0.920p 2 - 3.115,9 + 1 



(44) 



2. Strong coupling expansion 

In this case we identify (3 = g and a set corresponding to N u = 6 and N4 = 3, obtaining the solution 

q (-0.0000545(? 6 - 0.000489(? 5 - 0.00349g 4 - 0.0324 e 3 + 0.636p 2 - 1.555e + l) 
13 = ™ -0.327g 3 + L50V -2.180^+1 (45) 

In FigE] we have plotted P = —-jklnZ as a function of /3, as done also in Fig. 3 of [3l|. Our results show that 
both the Weak and Strong coupling expansions, resummed through the PPT converge to the exact result: this is 
particularly remarkable in the case of the weak coupling expansion. 

F. <f) 4 field theory in 1 + 1 

In a recent paper Nishiyama [32j has studied a lattice </> 4 model in 1 + 1 dimensions, described by the hamiltonian 



2** 



31 2 



(46) 



where i is the site index and 7^ and (pi are canonically conjugated operators. Notice that we have changed the notation 
in (32l | adopting the conventions used in this paper. 
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FIG. 9: (color online) Difference between the perturbative polynomial of order g 11 given by Nishiyama and the polynomial of 
order g 11 constructed with PPT working to order N u + Nd = 5. 



TABLE III: Comparison between the perturbative coefficients of [32|] and those predicted with PPT working with the set [3, 2]. 





br 


bs 


&9 


bio 


bn 


L exact 
On 

& [3,2] 

error (%) 


0.011061391245982 
0.011061133480144 
0.00233 


-0.0087493465269972 
-0.0087483769472128 
0.011081739435 


0.007096747591805 
0.007094602847397 
0.03022 


-0.005871428 
-0.00586767 
0.06403 


0.00493622 
0.00493037 
0.11851 



Using a linked cluster expansion Nishiyama has obtained the perturbation series in g up to order 11: 



E{g) = 0.66798625915577710827096201688+ 0.43100635014259473006095738275<? 

- 0.10148809521111863294125944502/ + 0.04803845646443637442034775341/ 

- 0.029018513979643624653232757064/ + 0.019777791330895673863274529570/ 

- 0.014454753622894705466341917665/ + 0.01106139124598227911409431586/ 

- 0.0087493465269972/ + 0.007096747591805/ - 0.005871428/° + 0.00493622/ 1 



(47) 



Since the perturbative series has a radius of convergence go w 1, an Aitken S 2 process was used in [32j to accelerate 
the convergence of this series. The accelerated series was then compared with the numerical results obtained using 
the Density Matrix Renormalization Group (DMRG), showing that the region of convergence could be enlarged up 
to g w 2. 

We will now apply our method to this problem and consider 



g{o) = e 



1 + 



(48) 



where as usual N u and Nd fix the asymptotic behavior for g — ► oo. As we do not know this behavior exactly we will 
work at order N u + Nd = 5 and take into account all the possible combinations of N u and Nd keeping the sum fixed. 
In this case only the coefficients b n with n going from to 6 are used, the remaining being a prediction of our method. 
In Fig|9]we have plotted the difference between the perturbative polynomial of order g 11 given by Nishiyama and the 
polynomial of order g 11 constructed with PPT working to order N u + Nd = 5. As one can see, the set [3, 2] provides 
the smallest difference, a result which suggests the asymptotic behavior of the energy as E oc yfg for g — > oo. 

In Table IIIII we have compared the exact coefficients calculated by Nishiyama with those predicted by the PPT 



,[3,2] 



/K x 



1 



from this results 



working with the set [3, 2]. The last row of this table shows the errors S„ = 100 x 

we can conclude that resummation through PPT allows to achieve a truly remarkable precision, the largest error 
being of about 0.1%. 
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FIG. 10: (color online) Comparison between the resummed energies to orders [3,2], [4,3] and [5,4] and the perturbative 
polynomials. 



In Fig llOl we have compared the perturbative polynomials for the energy from orders g 2 to g with the energy 
resummed with the sets [3, 2], [4, 3] and [5, 4]. There are several striking aspects which should impress the reader: first 
of all, the difference bewteen the three sets is extremely thiny, thus signaling that the convergence is extremely strong; 
in second place, the resummed energy confirms the DMRG result displayed in Fig. 2 of [321 ]: finally, the resummed 
energy is a multivalued function, with a branch point at g ps —1.025. This last finding is extremely interesting, 
because in [33] it was speculated the existence of a phase transition at g k, —2 (in our notation): the resummed 
energy plotted in Fig. 2 of [32j appears to have a singularity around g = — 1 (in our notation), i.e. in the same region 
where we observe the branch point 4 . Because of the use of a parameter g, PPT can deal with multivalued functions in 
a way which is not possible in conventional perturbation theory. Finally, the thiny dashed line in the plot corresponds 
to the numerical result obtained in [33[ using the Variational Sine Collocation Method (VSCM) within a mean field 
approach. 



G. Elastic theory 

Another example of divergent series has been studied in [34[. The authors of that paper have found out that the 
series for the inverse bulk modulus if as a function of the compression: 

— = b + hP + b 2 P 2 + ... (49) 
K 

has zero radius of convergence and they obtained an explicit expression for the coefficients in a simplified calculation. 
In the following we will uniform the notation to the conventions used in this paper and refer to the pressure P as the 
coupling g. The coefficients of the series are [34| 

& n = -(n + 2)%ta (50) 



where 



n/2 



/.-(-!,- r[^] ^ (61) 



4 Clearly, the branch point of function y = f(x) at a point x = xq manifests itself as a singularity in that point when it is calculated using 
the Taylor series around a different point. 
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TABLE IV: Comparison of the exact perturbative coefficient &7 for the elastic series with the approximate coefficient predicted 
by PPT using N u + N d = 5. b e 7 xact = - 7 -§^ « -27.2324646250573949498. The parameters are chosen /3 = A = a = y = l 
and a = 1/2. 





6^ 01 


b^ 


6 7 3 ' 21 


b^ 


b^ 


bf M 


error (%) 


-48.40832399 
77.75 


-27.41931538 
0.69 


-27.18605118 
0.17 


-27.21965479 
0.047 


-27.17382665 
0.21 


-27.46859099 
0.87 



TABLE V: Comparison between the perturbative coefficients of [34j and those predicted with PPT working with the set [2, 3]. 
We use ft = X — a = Y = 1 and a — 1/2. Coefficients with a dagger are input of the method. 





67 


bs 


69 


610 


fen 


rexact 
n 


- 27.2324646250574 


51.2814820354647 


- 100.257786099872 


203.060091610056 


-425.208317773323 


b [%3] 

error (%) 


- 27.2196547938467 
0.047 


51.1255824090448 
0.304 


- 99.224725761445 
1.03 


197.961957448658 
2.51 


-561.00565159656 
31.9 


error (%) 


- 27.2324646250574+ 



51.2814820354647 1 " 



-100.216793832003 
0.041 


202.536703574420 
0.26 


-579.771361475677 
36.35 



Here a is the surface tension, Y is the Young's modulus, a is the Poisson ratio, ft = l/fc^T and A is the ultraviolet 
cutoff of the theory. To apply our method we introduce 

9(Q) = 9 1 + E r lC " g " (52) 

and fix N u + Nd = 5. Working to this order the first predicted coefficient is 67. We have performed a calculation 
using ft = A = a = Y = 1 and a = 1/2. 



As we have seen from Table HVl the optimal set for N u + Nd is [2,3], corresponding to an asymptotic behavior 
1/K oc g°. At this order we have found 5 : 

q (-0.4923950887£> 2 - 0.94207053700^+ l) 
9 ^ 6 ' ~ 0.6331671415g 3 + 0.9634343699g 2 - 2.4724637506(? + 1 ^ ' 

Working with to order N u + Nd = 7 we have found that the best set is the [5,2], which provides a different 
asymptotic behavior for the inverse compression modulus, 1/K oc g 1 ^ 4 . In Tabic fVl we compare the predictions for 
the perturbative coefficients obtained with the two different sets: notice that the second set does not predict the 
coefficients 67 and b s . The second set gives more precise results than the first set for bg and 6 10 , but a slightly worse 
error for fen. 

In Fig llll wc have compared the perturbative polynomials of order 8 through 10 with PPT results corresponding 
to the sets [2, 3], [5, 2] and [6,6]: just as in the case of the lattice <fi 4 previously discussed wc observe a branch point 

in the resummed solution, corresponding to 3q 2 ' 3 ' ~ —0.412879, 3q 5 ' 2 ' ~ —0.401549 and ffg 6 ' 6 ' ~ —0.330171 with the 
different sets. 

If wc go back to the example of cf) 4 in zero dimensions, there we have seen that g = is a point where the function is 
not analytical and therefore the perturbative series is divergent. In the present example we can use the words of the 
authors of [34[ and say that "under stretching (g < 0) the true ground state is fractured into pieces. As a result g = 
cannot be a point of analyticity for K{g) and thus the series has zero radius of convergence". Our results however 
display a branch point not exactly at g = or close to it as in the case of the (f> 4 model in zero dimensions. 



Although the coefficients c n and d n are calculated exactly, we prefer to write them numerically to allow a more compact expression. 
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FIG. 11: (color online) Comparison between the perturbative polynomials of order 8 through 10 and the results obtained using 
the sets [2, 3] and [5, 2] with PPT . 



H. Elliptic integral of the first kind 



Consider the elliptic integral of the first kind: 

r^ 2 dt 



r i 

E(g) = K(g) = / -= 
Jo \/l 



yl - g sin 2 1 

which diverges for g — > 1. It also obeys the series representation 



(54) 



^) = fE%^V (55) 

fc=0 

which converges for \g\ < 1. 

We want to show that it is possible to resum the perturbative series using the PPT. We choose the functional form: 

9iQ) = 9 1+ St Cngn ■ ^ 

The choice of N u + 1 = Nd is not arbitrary: with this choice we have that lim^oo g(g) = g < oo, which means that 
the resummed function will have a singularity precisely at g = g. 
Using N = 2 we find 

/ 38 lg 2 , 187£ , i 

f l SRS«in ~ 9940 ' 1 



K 35840 1 2240 

9{8> ~ 2301g 3 , 1181g 2 , 1447g , ] ' ' 

286720 8960 " + " 2240 + 

which predicts the singularity of the elliptic integral at 

- g w = ™> w L324 

767 

Increasing iV this singularity moves towards its exact value, g = 1; for example, using JV = 3,4 and 5 and find 

g™ = ^ „ f.332 , 5^ 5 1 = *™ « 0.670 , 5 M = 69944792 „ 1.035 . (59) 
y 3752 ' y 8216 ' y 67596985 v ' 

The reader will notice that the singularity predicted by the set [4, 5] falls below the exact singularity: the reason 
for this behavior is easily understood looking at Fig|T2] As a matter of fact the set [4, 5] (the thin line in the plot) 
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FIG. 12: (color online) Comparison between elliptic integral K(g) and the PPT approximations obtained with sets [2, 3], [3, 4], 
[4,5] and [5,6]. 



has a branch point close to g — 1, and therefore the singularity belongs to the nonphysical branch. Notice that the 
set [5,6] provides an excellent approximation. 

Let us now compare the expansion of K(g)^ 2 ' 3 ^ around g = with the exact result, provided by the series (|5"5")) . We 
have 



K{g) 



[2,3] 



and 



K(g) 



7T 



7T ng 9tt/ 25tt/ 1225tt/ 3969tt/ 53361tt/ 206126367tt/ 
2 + ~8 + 128 + 512 + 32768 + 131072 + 2097152 + 9395240960 
405813405891tt/ 810831328918663tt/ 1639189758117069059tt/ 

21045339750400 + 47141561040896000 + 105597096731607040000 
33455928294943808886877rg n 6882636481373124653844491tt/ 2 

236537496678799769600000 + 529843992560511483904000000 + "' 
1.5707963267949 +0.392699081698724 5 + 0.220893233455532/ + 0.153398078788564/ 
0.1 17445404072494/ + 0.09513077729872/ + 0.079936278146842/ + 0.06892479746239/ 
0.060578751866013/ + 0.054035158997418/ + 0.0487671220263656/ 

0.0444347725101475 5 n + 0.0408090692936214/ 2 + . . . (60) 
irg 9tt/ 25tt/ 1225tt/ 39697r ff 5 53361tt/ 184041tt 5 7 



2 8 128 512 32768 131072 2097152 8388608 
41409225tt/ 147744025tt/ 213342372l7r.g 10 777553604l7rg n 45702872952l7r 5 12 
+ 2147483648 + 8589934592 + 137438953472 + 549755813888 + 35184372088832 + "' 
w 1.5707963267949+0.392699081698724.9+0.220893233455532/ + 0.153398078788564/ 
+ 0.117445404072494/ + 0.0951307772987205/ + 0.0799362781468415/ + 0.0689246479939603/ 
+ 0.0605783039009417/ + 0.0540343513190498/ + 0.0487660020654424/° 

+ 0.0444334853530168 5 n + 0.0408078363745588/ 2 +.. . . (61) 

Notice that coefficients starting from 7 and higher, are predictions of the PPT: we see, for example, that the 
coefficient of the term of order g 12 is predicted with an error 0.003%! 

I. Virial coefficients 

As a last example we consider the low density virial expansion of the pressure 

P 



k B T 

n— 1 



E B ^ fc ' ( 62 ) 
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TABLE VI: Virial coefficients for a hard spheres in 2 and 3 dimensions given in Table I of [35| and predictions using PPT with 
different sets. 















D -■ 


= 2 


D -■ 


= 3 


Ref.f35l 


0.0362193 


0.0199537 


0.0013094 


0.0004035 


[0,5] 


0.03739998 


0.02496595 


0.0023400 


0.0031580 


[1,4] 


0.03625994 


0.02008503 


0.0013509 


0.0004884 


[2,3] 


0.0362321 


0.0199843 


0.0013165 


0.0004198 


[3,2] 


0.0362551 


0.02006717 


0.0013404 


0.0004664 


[4,1] 


0.0368599 


0.02258546 


0.0017325 


0.0014442 


[5,0] 


0.1747048 


0.75984885 


0.0222648 


0.0735264 



where the Bk are the virial coefficients and p is the density (not to be confused with g). Table 1 of [331 contains the 
numerical values of the virial coefficients for hard spheres in D dimensions, with 2 < D < 8. In the following we will 
uniform the notation in (|62|) to the notation adopted in this paper and call g the density. 
As usual we adopt the functional form 

9{q) = 9 1+ S: lCngn ■ ^ 

In Table IVII we have used PPT with N u + Nd = 5 to predict the virial coefficients Bg and Bio from the previous 
one. As one can see the set [2, 3] provides highly precise results. We have then used the set [3, 4] which has the same 
behaviour for g — > oo to obtain an estimate for the eleventh virial coefficient. To orders [2, 3] and [3, 4] we have found 

[Bii/fla ] 18 ' 31 = 0.01094432 , [B u /B 2 10 ] [M = 0.01090061 . (64) 

Since these results are not (yet) available in the literature, the prediction made here will be a strong test of the 
present method once the calculation of B\\ will be made. In Tabic IVIII wc compare our predictions for the virial 
coefficients going from B\\ to Bis with the predictions made in [35]. Our predictions are very close to those made by 
Clisby and McCoy for D = 2. 

Notice that finding N u + 1 = Nd has an important effect: as discussed in the previous example of the elliptic 
integral, in this case the solution will have a singularity at a finite value of g. For example, if we consider the set [3, 4] 
the tranformation reads 



q (0.3004g 3 + 1.8876p 2 + 2.71410 + l) 
0.2584g 4 + 2.1266e 3 + 4.6986g 2 + 3.7831g + 1 



9(0) = n n e o„„4,, n 1 .3 I a .no.J I o ,0,1 I I 1 ( 65 ) 



and the singularity is predicted to be fall at 

g [3A] = lim g(g) w 1.1625. (66) 

Q — >00 

We would like to stress that this singularity is "physical" , i.e. it is a singularity of the resummed function, in 
contrast with the singularity falling at the radius of convergence of a series. We can use the previous example of the 
elliptic integral to better understand this point: in that case the perturbative series around g = had a radius of 
convergence 1, coinciding with the location of the true singularity of K(g). 

If we trust our result, we may conclude that the becomes infinite at a finite density g ~ g^' 4 \ 
We have also considered the virial series in D = 3 dimensions. Also in this case we have found out that, working 
with N u + Nd = 5 the optimal set corresponds to [2,3] (see Table [VT|) and therefore the virial series is expected to 
have a singularity at finite density: 

g\m = i im g ( Q ) ~ 1.43439 . (67) 

Q — >00 

In Table Ivnl we have also compared our predictions obtained with the set [3, 4] for D = 3 with those made in [35| . 
Unlike in the previous case, our result agree to some extent with those of [351 ] only for the coefficient i?n, whereas 
completely different predictions are made for the remaining coefficients. 
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IV. CONCLUSIONS 

We have developed a new method, Parametric Perturbation Theory (PPT), which is alternative to the ordinary 
perturbation theory, i.e. does not amount to an expansion in any physical parameter. Wc have shown that PPT can 
used either as a fully autonomous perturbation scheme, as done in Section [HI or it can be applied to the coefficients 
of the perturbativc expansion, resumming the series and providing physically meaningful results, as done in Section 
IIIII There are several aspects of our method which should make it very appealing: 

• since PPT can use perturbative results as an input, it can be applied with limited effort to the huge amount of 
problems which have been studied perturbatively; 

• it provides analytical approximations; 

• unlike variational methods, such as the LDE or VPT, our method does not require any optimization in a 
variational parameter; 

• although the asymptotic (strong coupling) behavior of the solution can be used, when known, to refine the 
functional relation g = g(g), PPT is capable of selecting the most appropriate asymptotic behavior of the 
solution within a class of different behaviors allowed to a given order; 

• it predicts the unknown perturbative coefficients with high precision; 

• it can easily describe multivalued functions and therefore is capable to produce branch points at finite order, as 
observed in the examples: if these points are related to phase transitions of a system, as claimed in (32lj in the 
case of <f) A in 1 + 1, then our method could provide an alternative tool to the study of critical phenomena; 

• it can produce singularities in an observable working at finite order, as seen for the cases of the elliptic integral 
of first kind and for the virial coefficients of a hard sphere gas; 

• most importantly, it can produce the nonperturbative imaginary part of an observable, which appears when a 
system becomes metastable. 

Future directions of work will certainly include the development of an autonomous perturbation scheme for quantum 
mechanical problems, in analogy to the one developed in classical mechanics and the application of our results to 
resum perturbative calculations in quantum field theory. It will be also interesting to apply PPT to obtain new 
analytical approximations for special functions of relevance in Physics, as done in this paper with the elliptic integral 
of the first kind. 
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TABLE VII: Predicted coefficients for approximants with 10 exact coefficients for D — 2 and D = 3. Comparison between the 
predictions of [35| and the predictions obtained using the set [3,4]. 





Bn/Bl° 


B12/S2 1 


B13/B 42 


B14/B}, 3 


B15/B 14 


S 16 /B 2 15 




Bis/Bl 7 


D = 2 


















Ref. [351 


1.089 x 10~ 2 


5.90 x 10~ 3 


3.18 x 10~ 3 


1.70 x 10~ 3 


9.10 x 10~ 4 


4.84 x 10~ 4 


2.56 x 10~ 4 


1.36 x 10~ 4 


[3,4] 


1.0901 x 10~ 2 


5.9235 x 10~ 3 


3.2117 x 10" 3 


1.7421 x 10" 3 


9.4698 x 10~ 4 


5.1638 x 10~ 4 


2.8247 x 10" 4 


1.5492 x 10" 4 


D = 3 
Ref. [351 

[3,4] 


1.22 x 10~ 4 
1.1599 x 10~ 4 


3.64 x 10~ 5 
2.2229 x 10~ 5 


1.08 x 10~ 5 
-8.5616 x 10~ 6 


3.2 x 10~ 6 
-1.8088 x 10~ 5 


9.2 x 10~ 7 
-2.0325 x 10~ 5 


2.6 x 10~ 7 
-2.0112 x 10~ 5 


-1.9136 x 10~ 5 


-1.7971 x 10~ 5 



